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j^ I We study with lattice Monte Carlo simulations the interactions and macroscopic behaviour of a large 

number of vortices in the 3-dimensional U(l) gauge+Higgs field theory, in an external magnetic 

Cl^l field. We determine non-perturbatively the (attractive or repelling) interaction energy between two 

I5 ■ or more vortices, as well as the critical field strength He, the thermodynamical discontinuities, and 

the surface tension related to the boundary between the Meissner phase and the Coulomb phase in 
the type I region. We also investigate the emergence of vortex lattice and vortex liquid phases in the 
type II region. For the type I region the results obtained are in qualitative agreement with mean field 
theory, except for small values of H^, while in the type II region there are significant discrepancies. 

C^ I These findings are relevant for superconductors and some models of cosmic strings, as well as for the 

electroweak phase transition in a magnetic field. 
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1 Introduction 

The continuum 3-dimensional (3d) U(l) gauge+Higgs field theory (scalar electrodynamics) is of in- 
terest for several reasons. First of all, it is the Ginzburg-Landau (GL) theory of superconductivity, 
believed to be applicable for low- Tc and possibly also for some aspects of high- Tc superconductors. 
Second, it is an often-used toy model for studying line-like cosmological topological defects, strings 
(adopting here the condensed matter terminology, we will mostly use the word vortex for these de- 
fects) . Finally, the existence of topological excitations makes the model interesting also as a theoretical 
playground: much of the terminology related to our understanding of confinement, for instance, arises 
from superconductors. It is rather remarkable that for this theory, the topological observables of the 
continuum limit can be easily generalized in a gauge invariant way to the case with a finite lattice 
cutoff 1,1. 

In all these different contexts, it is of interest to ask how the behaviour of the system changes when 
an external magnetic field is added. In superconductors, an external magnetic field can be imposed 
experimentally, and it leads to the emergence of completely new phases, displaying for instance broken 
translational invariance. Due to magnetohydrodynamic diffusion, essentially homogeneous magnetic 
fields could also exist in cosmology and be relevant for the phase transitions appearing there; in any 
case, the core of each long U(l) vortex carries magnetic flux. From the theoretical point of view, a 
small external magnetic flux allows also to probe the topological properties of the zero-flux system in 
a gauge-invariant and systematic way. 

In a previous paper |3| , we showed how one unit of external magnetic flux can be imposed gauge- 
invariantly on the lattice, and how it can be used to determine non-perturbatively the free energy per 
unit length, i.e. tension, of a long vortex. The vortex tension thus defined was shown to constitute a 
gauge-invariant (but non-local) order parameter for the system, differentiating between the symmetric 
( "Coulomb" , "normal" ) and broken ( "Meissner" , "Higgs" , "superconducting" ) phases. In the present 
paper, we study the case of a larger magnetic field. We show how to measure the free energy of a 
system of two or more vortices, thus seeing whether vortices attract or repel each other. This allows to 
divide the phase diagram into the region of type I and type II superconductors. We also increase the 
magnetic field further on and inspect whether this can lead to the emergence of new phases, triggered 
by the interactions of vortices, which are themselves macroscopic objects from the point of view of 
the original continuum quantum field theory. 

In condensed matter physics vortices, of course, are a thoroughly studied topic (see, e.g., |Q-0])- 
The emphasis is somewhat different, though. There the existence of a multitude of different scales 
and parameters implies that even the starting point, theory or model, is not uniquely known. The 
multitude of scales is related to the fact that one usually stays away from the immediate vicinity of 
the transition point, as a result of which the mean field approximation also often works very well. In 
particular, if one is studying the GL theory, one often approximates it by a simpler theory (the frozen 
gauge model, the XY model, etc.), using the narrowness of the fiuctuation region as an argument. 
These arguments may be correct for physics, as it is only experiment which decides which effective 
theory is best, but they are not correct for the GL theory as such. 

In the particle physics context, in contrast, given definite symmetry principles and the requirement 
of ultraviolet insensitivity, the form of the theory is uniquely defined; in the case of the U(l)-|-Higgs 
theory, keeping only relevant operators, it depends just on two dimensionless parameters. Moreover, 
close enough to the transition point, mean field approximation completely fails, and the theory is a 
genuine quantum field theory. It is quite striking that even for a simple theory such as U(l)-|-IIiggs, the 
properties of the phase diagram are not yet completely known. We hope that studying the response 
of the system to an external magnetic field will shed new light on this issue. 

Somewhat similar topics (with different methods and in 4d) were first studied with lattice simula- 
tions in g]. External magnetic fields in the 3d SU(2)xU(l) theory (corresponding to the Standard 
Model) were studied in Q, and in the 3d U(l) theory with fermions instead of scalars, in ||lO| . 



The plan of the paper is the foUowing. In Sec. ^ we define the theory in the continuum and on the 
lattice, and review how an external magnetic field can be imposed and how the free energy related 
to vortices can be measured. In Sec. ^ we review the mean field results for the vortex and surface 
tensions and for the behavior of the system in an external magnetic field. In Sec. Hwe consider the 
full path integral and discuss the validity of the mean field approximation. In Sec.^ we present our 
simulation results for the type I region, and in Sec. |g for the type II region. We conclude in Sec. 0. 



2 An external magnetic field and vortices on the lattice 

The theory in continuum without external field. The theory underlying our considerations is 
the U(l)+Higgs theory, corresponding to high temperature scalar electrodynamics on one hand |ll|| - 
flil and to the Ginzburg-Landau model of superconductivity on the other jlj] . It is defined by the 
functional integral 



Z = /'l?A,I?</.exp[-5(A„ (/.)], 
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djAi and Di 



proper powers of the scale Cg, the theory can be parametrized by the two dimensionless ratios 



di + ie^Ai. When all dimensionful quantities are expressed in 
be p 

.= ^, . = ^, (2.3) 

where m\{^) is the mass parameter in the MS dimensional regularization scheme in 3 — 2e dimensions. 
For completeness, one may note that the standard dimensionful textbook coefficients a,b [|l^ of the 
Ginzburg-Landau free energy and the dimensionless GL parameter k are related to y, x by 

mc^ 1 /mr\ 2 



y 



where a is the fine structure constant and m is an effective mass parameter. Note the huge numerical 
ratios entering here. 

The theory on the lattice without external field. The discretized lattice action corresponding 

PgY. J4(x)- ^^Re0*(x)f/,(x)0(x + z) 

(2.5) 



to Eq. (2.2) is 
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where ai(x) — ae3Ai(x), aij(x) — ai(x) + aj(x + i) — Q:i(x + j) — aj(x), C/i(x) = exp[iai(x)], and 
where the lattice couplings /3gj P2 are related to the continuum parameters and the lattice constant a 
by um (for C'(a)-corrections, see |l7|) 
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Note that we must use here the non-compact formulation for the U(l) gauge field, to avoid topological 
artifacts (monopoles) at finite a (with the compact action one needs to go closer to the continuum 
limit to start with than /3g > 1 as we have here, at least to (ic > 4; see [ll2[). 

The thermodynamical ensembles v^rith an external field. What we do in this paper is an 
extension to n vortices of what was done for one vortex in B. Physically, we impose an external 
magnetic flux $b through the volume, 

e3$s = eaSLiLa = Svrn. (2.8) 

By convention, we have fixed the flux density as B = (0,0,5), and use an N1N2N3 lattice with 
Li = NiU, V = L1L2L3, A — LiL2- On the average, after taking into account translational invariance, 
this flux density determines that {F12) = B, although in individual configurations the flux may not 



be evenly distributed. Together with the requirement of periodicity of observable quantities, Eq. (2.8) 
guarantees that there are n vortices passing through the system, and we want to measure the free 
energy F{B) of such a conflguration. Technically, flxing the flux is equivalent to multiplying the 
integrand in Eq. (|^) by the delta-function 

n^l^^sU dxidx2e:iFi2{^) - 2^n) , (2.9) 

where the product has been written in a discretized form. 

We shall use interchangeably the notation F{n) or F{B) for the free energy thus obtained, noting 
that n and B are related by Eq. ( |2.8| ) and that F also depends on the parameters e§, j/, x (see Eq. (|2.3D) 
of the theory. Note that in our convention the free energy is dimensionless, Z = exp(— F), so that it 
corresponds to the free energy divided by temperature in standard terminology. 

In a thermodynamic sense, flxing the flux means choosing a microcanonical ensemble. The results 
can also be transformed to the canonical ensemble G{H), 

G{H) = F{B) - VHB, F'{B) = VH, (2.10) 

where V is the volume of the system and H is the external held strength. Technically, the canonical 
ensemble means multiplying the integrand in Eq. (EHJ) by the factor 



exp[+H f d^xFui^)]. (2.11) 

Imposing the magnetic fiux on the lattice. To impose a flux ^b of n units of 27r/e3 on the 

lattice, we choose here to proceed as follows. We deflne a modifled theory by 



^('") = / T)r^T)rh p-S[a,0;m] 



(2.12) 

where the degrees of freedom are the periodic link angles 

— 00 < ai{x.) < 00, ai{Ni + 1,X2,X3) = ai{l,X2,X3), etc., (2-13) 

the periodic scalars (/)(x), and 

X 

= S[a,<j>]+l3G^{2Trmai2{xo,yo,X3) + 27T^m^). (2.14) 



Since, to leading order in a ^ 0, 

Oii2{xo,yo,X3) = a'^e3Fi2{xo,yo,X3) 



63^1 



(2.15) 



we have forced a flux of — 27rTO/e3 through the lattice; the stack of plaquettes parallel to the X3 axis 
at the position {xo,yo) is, in fact, a "Dirac string" c arryi ng the flux 27rTO/e3 in the —X3 direction. 
A crucial fact now is that, due to the periodicity ( 2.13 ), 



y^ ai2{xi,X2,X3) = 0; 



any 0:3, 



(2.16) 



and thus the total flux through any plane vanishes identically. This implies that the flux +27Tm/e3 
must return through the system in the +X3 direction but now in a manner specified by the dynamics 
of the theory. This response is the object of the study here. 

The Dirac string has been imposed by modifying one special plaquette on each 0:3 plane leaving the 
part of the action involving scalars unchanged. When m is an integer, integrating over all periodic 
field configurations makes the path integral ( 2.12| ) translation invariant. One could also, equivalently, 
impose the flux by making a non-periodic change in the boundary conditions for the link variable a^, as 
was done in Ref. H|. Although translation invariance is broken with non-integer m, the path integral 



(2.12) is well-defined, and non-integer values will indeed be used to interpolate between integers. 

Finally, in principle one might also attempt to perform simulations with the canonical ensemble, as 
was do ne in |q| . In our case this can be accomplished by coupling an external field to the flux according 
to Eq. ( 2.11 ), i.e., adding a term HL3 J (PxFi2{x) = HL32im/e3 to the action. This would promote 
n to a dynamical variable, for which only integer values are allowed. However, in our case it would be 
very difficult to obtain an efficient update for n. This is because the update is of semi-global nature: 
by performing it one attempts to make a large change in a configuration without any interpolation. 



Measuring the free energy. The aim now is to measure the change in the free energy caused by 
switching on the magnetic field, F{n) — F{n = 0). Since one cannot measure the absolute value of 
the free energy on the lattice, we eliminate the unknown i^(0) by taking a derivative with respect to 
m and integrating back. This leads to the following result for the free energy per unit length (always 
relative to F(0)): 



Lsej 



= 2tt^(31 



dm W{m) 
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(2.17) 



The subscript m emphasises the fact that the expectation value is computed using the action (2.14) 
with the defect. We recall that we have chosen to express every dimensionful quantity in units of an 
appropriate power of 63. For instance, Fin)/ L3 was made dimensionless in Eq. (2.17) by dividing it 
byei- 

Since W{m) is the quantity measured in our simulations, it is useful to have a feeling of its mag- 
nitude. It is essen tially the expectation value of the stack of plaquettes at the position of the Dirac 
string. Eq. ( 2.14 ) for the first implies that the action is small for 012(2^0, 2/0, 2:3) = —2T:m. This 
large negative contribution cancels the term 2m in ( 2.17 ). This rough estimate can be improved by 
evaluating W{m) including only the gauge field in the action and using ( 2.16 ). The plaquettes ai2 
are of two types: one at the string and N1N2 — 1 not at the string. Using shift symmetries of the 
action one can then prove that 



Type I 



Type II 
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Figure 1: The behaviour of the intersection of flux hues with increasing total flux on a flnite lattice in 
type I and II cases (see text). A dark area is a region of the symmetric phase carrying the flux, while 
a white area is a region of the broken phase without a flux. The notation is somewhat symbolic: in 
the type I case, the flux lines are in reality on top of each other. 
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B given by (|2.8|). In fact, the shift symmetries 



for all 771, which simply corresponds to F 

also hold for the full theory if ni ^ n = integer, so that Eq. (2.18|) will give W{n) exactly when n 



integer. The behaviour of W{m) between integer values is seen from the numerical computations in 
Figs. |(a-c). 

Measuring the field strength. In addition to the total free energy of n flux units, we will be 
interested in the increase of free energy when the flux is increased by one unit, dF(n) = F{n)—F(n—1). 
Using dB ■ L1L2 = ^Tr/e^ and Eq. ( 2.10| ), this is, in fact, a measurement of H: 
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3 Mean field results and 

the determination of physical observables 

In this section, we discuss the behaviour of the system at mean field level, and how this pattern can 
be employed in determining some of the physical observables of the full quantum theory. 

One or two vortices. In the mean field approximation the phase of the system is determined by 



the configuration that minimizes the action (2.2) under the boundary conditions imposed. In the bulk, 
for zero flux the system has two phases depending on y; a broken phase at y < and a symmetric 
phase at y > 0. The transition between the phases is of second order. In the broken phase the scalar 
and vector masses are 



= y^ei, mr = ^^el (3.1) 

Assume now that the whole system is in the broken phase and ask what happens if one starts 
increasing the magnetic flux in the 3:3 direction keeping the transverse area fixed (Fig. |^). This is 
concretely how our simulations will be organised. 

Taking first n = 1 one vortex line appears somewhere. By definition, the vortex tension T is the 
free energy of this configuration divided by the length of the system in the xa-direction. In the mean 



field approximation, T can be calculated by minimizing the action (2.2) with cylindrically symmetric 
boundary conditions. The result is |lq 

^(^^= 1) _ .2 f-y\ _c/./^^ _T.MF 



T = -^ ^ ^ei{^] tt£{V2x) = T^^^ , (3.2) 

L3 \ X J 

where £{l) = 1, £(0.3043) = 0.63 {x = 0.0463) and £(2) = 1.32 {x ^ 2). 

When n is increased to n = 2, a difference between type I (a; < 1/2) and type II {x > 1/2) appears: 
in the former case vortex lines attract and in the latter they repel each other. If we form the quantity 

Fin = 2) , , 

T. = ^^, (3.3) 

we thus expect that for type I, 

r-T2>0 (type I), (3.4) 

due to the binding energy of two vortices. This also holds in the thermodynamic limit. For type II, 

T-T2<0 (type II), (3.5) 

where the inequality holds on finite lattices and the equality in the thermodynamic limit. Thus T — T2 
is an order parameter for separating the type I and type II regions. In ||j it has been shown that 
T itself, on the other hand, serves as an order parameter for the normal <-> superconducting phase 
transition in this theory: it vanishes in the thermodynamic limit for y > ydx). 

The boundary between type I and II, x = 1/2, is mathematically very interesting since the classical 
equations of motion then simplify considerably. The dynamics of n vortices can then be largely 
solved ||l^. We shall not study this special case here. 

3.1 Type I, x < 1/2 

The general pattern. Consider then type I and increase n further. The lines parallel to X3 first 
cluster on the a;i,a;2 plane as a 2d bubble (= 3d cylinder) of symmetric phase, but at some stage it 



becomes more favourable to form a slab (see Fig. |l|). The cylinder and slab have interfaces separating 
the broken and symmetric phases. The tension of this interface can be written as 

_ MF (^?-//2) / \ 4 /Q a\ 

a — a = s(a;je3, (j-d) 

where s{x) is the extremal value of an integral already written down by Ginzburg and Landau. For 
small X 1221 , 

s(a;) = ^-1.24x1/4 + ... (3.7) 

and one can analytically prove that s(x) vanishes for x — 1/2 and becomes negative for x > 1/2. For 
intermediate values numerical methods have to be used pQ, El| . For the value of x used in the type I 
simulations we find s{x = 0.0463) = 0.643. 

With continuously increasing n the slab gets thicker. The broken phase region gets smaller and, 
ultimately, it reduces to a 3d cylinder of the broken phase, a hole in the magnetic field configuration. 
At some critical value of the flux, corresponding to a critical value of the field strength He, the system 
goes entirely to the symmetric phase. 

The detailed finite volume behaviour. Consider now an inhomogeneous configuration in which 
the area As is in the symmetric phase, and the area A — As is in the broken phase with no B. Let us 
denote the broken phase action density as 

~™^^-V^ He^nr^^el (3.8) 



4A3 " 2 ^' ^- ' V2^ 

In the symmetric phase, the magnetic flux density is Bs = BA/Ag, and the action density is ■^B^. If 
dAs is the length of the boundary, the action per length is 

£ = ^AsB's + iA- As){-^H'J + adAs, (3.9) 

where the three stages discussed above correspond to 

cylinder of symmetric phase: Ag ~ irr^ , dAg = 2Trr, 

slab: As = Larf, dAs = 2^2, (3.10) 

cylinder of broken phase: As — A — nr'^, dAg = 2Trr. 

These now have to be minimised for cylinder radii r and the slab width d (it was assumed that 
L2 < Li) under the condition $ = BA = BsAs- The configuration corresponding to the minimum of 
F/L3 then is the stable one. The minimisation is simplest for the slab: then 

d^^L,. (3.11) 

For the cylinders one has to minimise numerically, in general. However, an analytic result is obtained 
if the surface term is so small that it can be neglected in the minimisation, though it affects the value 
of F/L3. This is true if ct ^ H^r at the minimum value of r. Then 



-Fcyl.symm. _ f rj rj ^ tt2\ t r , r, _ JTtBLiL2 



= iJ,i3--i7,^ iiL2 + 2aW^-^. (3.12) 



L3 \ 2 / V Hr 

For a slab configuration with L2 < ii, one obtains 



^ = (hcB - ]^Hl ] L,L2 + 2aL2. (3.13) 



Finally, for large n one has a cylinder of broken phase within symmetric phase and 

-^cyl.brok. 



1 



= [H,B--HnL,L2 + 2aJ7rLiL2[l- — 



B 



Comparing (3.12) and (3.13) one sees that a cylinder is favoured if 

B L2 _ Bi 
He nil He' 



(3.14) 



(3.15) 



To see when the transition to the homogeneous symmetric phase takes place, one should compare 
with its free energy: 



One sees, comparing ( 3.13 ) and ( 3.16 ), that this is smaller than even the slab free energy if 



B>Hr-2j— = B2. 

^1 



(3.17) 



The transition to the symmetric phase will take place directly from the slab stage if the free energy 
of the 2nd cylinder stage is at B2 larger than that of the slab: 



He Li 

2-K L2 V Li 



(3.18) 



Our simulations will, in fact, satisfy this condition. To see the stage with the cylindrical hole in the 
magnetic field configuration, one has to go to still larger lattices. This result is also confirmed by a 
full numerical minimization of Eq. (3.9). 

The determination of He and a. On the basis of the above equations we thus have a definite 
scenario for what happens when we increase the flux n further from n = 2 at finite volumes, and how 
this is related to the physical properties of the transition in the thermodynamical limit. 

To see what Eq. (3.12) predicts for the measured H{B) (Eq. (2.19)), simply take a derivative of 
Eq. (3.12) with respect to B (we recall that the relation of B and n is as in Eq. ( |2.8| )). Then, for 
small B the measurements should behave as 



HiB) =He + (T 



LiL2BHe 



(3.19) 



For B > HeL2/nLi — Bi (cf. Eq. ( 3.15| )), the slab configuration dominates and leads to a constant 
plateau, 

H{B) = He. (3.20) 

Finally, when B = B2, the system enters the Coulomb phase in which 

H{B)=B. (3.21) 

Thus, we can determine He from the plateau according to Eq. (3.20), and then a from Eq. ( 3.19| ). 



The determination of thermodynamical discontinuities. We have here considered a system 
with a fixed S, on which the magnetic field H depends, but it may be more natural to move to the 
canonical ensemble in which H is fixed and B is allowed to fluctuate. This means simply inverting the 
relation of B and H. From this point of view, we simply have a standard first-order phase transition 
at H — B2, and when H is below that, B vanishes. The configurations with a cylinder or a slab 
describe the system on the transition line. 

The thermodynamical discontinuities related to the first-order transition are directly given by the 
properties of the canonical free energy G{x, y, H/e^). The function G{x, y, H/e^) is continuous across 
the phase transition, but its derivatives are not: 






A^^^=etA{^*<j)), A^^^^-elAB. (3.22) 



For fixed x, the latent heat L of the transition is defined as the discontinuity in the "energy" variable E, 
obtained from G by a Legendre transformation with respect to y,H/e'^: 

L = AE = -yA— - H,A— = V [-ye\A{(^*^) + H,Ab\ . (3.23) 

Note that for fixed x, the identity AG'(x, y, He) — leads to the Clausius-Clapeyron equation, relating 
the different discontinuities: 

^a^— ^^ai? ^ etAir^)^—AB. (3.24) 

Thus it is enough to measure one of the discontinuities, and the curve Hc{y)- Below we choose to 
discuss AB (on finite volumes, AB = B2 instead of AB — He). Finally, let us recall that at the mean 
field level Eq. (p. 24) is satisfied through 



(^)el<0, H, = -^el AB = H,. (3.25) 



. 2x 

3.2 Type II, x > 1/2 

The determination of vortex tension and vortex interaction energy. In the type I region, 



the vortex tension can be determined directly from Eqs. (3.2), (3.3). In the type II region, the issue 
is more involved, since there are large finite volume effects. Indeed, as discussed after Eq. (3.5), the 
difference between T, T2 vanishes in the thermodynamical limit. However, this finite volume behaviour 
is well understood, and can be employed for the measurement of the properties of the vortex system. 

The basic idea is that due to the repelling interaction, the vortices tend to form a vortex lattice in 
the type II region. In fact, due to the periodic lattice boundary conditions, even one vortex forms a 
square lattice with its periodic counterparts. Thus one can use finite size scaling with n = 1 to study 
the interaction energy of two vortices. The method can be continued with n = 2, etc. Of course, 
this is not equivalent to a real vortex lattice, in which the distances between the vortices fluctuate; 
here the distance is fixed by the lattice size. When n grows to much larger values so that the vortex 
separation is no longer determined by the lattice size, the vortices are instead supposed to form a 
physical (triangular) Abrikosov lattice. 

Mean field theory offers an ansatz for the interaction energy: if d 3> l/my is the physical distance 
between two vortices (it can be the lattice size L in finite size scaling studies of T or T2 or the average 
distance between vortices, d ~ \J LxL^jn = ^2-k jiBe-^, at large n), then the interaction energy 
is 11 

£12 = 2TXKo{mvd), (3.26) 



where X = l/£{\/2x), Ko{z) is a Bessel function (~ ■y/7r/(2z) exp(— z) for z -^ cxj), and mv is the 
photon mass (in the broken phase). We have chosen to parametrize the prefactor in terms of the 
tension T for later notational convenience, although, in fact, the physics of the prefactor is different 
from that in T: the prefactor is only sensitive to asymptotic fields and the photon mass, and thus 
£{y/2x) cancels bet ween the mean field expressions for T and X . Let us also recall that the function 
Kq appears in Eq. ( 3.26 ) because it determines the profile of a single vortex. 

In the following, we will need an ansatz for the total free energy of a system of vortices forming 
either a square or a triangular lattice. We define the lattice sums 



Ssiz) 
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ni — — oc 
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St{z) = Y. Ko[z^J{2/V3){nl+nl + nin2)), 



(3.27) 



where the term rti = rt2 = is excluded and where the factor 2/\/3 follows by requiring that the two 
lattices have the same number of sites within a large rectangular area. (From the point of view of 
a triangular lattice, it would be more natural to use a lattice geometry which is not rectangular but 
has Li : L2 ~ 2 : v3, but we stick here to Li — L2.) A numerical computation of the sums shows 
that they are very close to each other in the relevant range of 2; this is the known fact that it is 
energetically difficult to distinguish between a square and a triangular lattice. This also means that 
it is quite difficult to differentiate between a lattice structure and a liquid phase, where the positions 
of the vortices fluctuate but have suitable average distances; we do not introduce a different ansatz 
for the liquid case. 

We now parametrize the ansatz for the free energy of n (n ^ 1) type II vortices in an approximately 
triangular configuration measured on a square lattice {L = Li = L2) as 






T 



n[l+XSTimvL/V^)] , 



(3.28) 



where we allow eventually T, X, my to deviate from their mean field value. For n = 1, we have instead 
one vortex interacting with its replicas on an infinite periodic square lattice of period L ^ Li ~ L2, 



T{L) 



F{n= I) _ T 
e 



^ 2 =-[l + XSs{mvL)]. 



-365 



(3.29) 



and this can be used for the finite size scaling studies of the tension T and the interaction energy X. 
For two vortices the period is L/^/2, 



T2{L) F{n ^2) _T 



-365 



1 + XSs{mvL/^/2) 



(3.30) 



which allows similarly to determine T,X. 



From Eq. ( 3.28| ) we can further, by differentiating with respect to n in analogy with Eq. (2.19) 
derive that 



H 

73 



evaluated at z = mvL/^/n. 



dF{n) 



T 

2^, 



1 + X[St{z)-IzS^{z: 



(3.31) 
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The determination of Hci,Hc2- When the volume and the number of vortices are large, the 
question is what is their equilibrium configuration. Let us note first that at small enough H, one is 
again in the Meissner phase just as in the type I r egion . The critical value Hd for the transition to 



the vortex lattice phase can be obtained from Eq. ( 3.31 ): in the infinite volume limit, St{z) — > and 



% = -^W-A5(V2^). (3.32) 

The magnetic flux B increases continuously from zero, unlike for type I where there is a discontinuity. 
Once H > Hd, the vortices form a triangular Abrikosov lattice structure on the mean field level. 
Increasing H further on, the lattice structure finally disappears when the energy squared corresponding 
to the lowest scalar Landau level becomes positive. This happens at 

H^^^-yel (3.33) 

Note that the physics determining Hd , Hc2 is quite different from that determining H^ in the type I 



region, Eq. (3.8); at a; ^ 1/2 the expressions of course meet. 



4 Fluctuations 

The mean field approximation discussed above misses some essential features of the full path integral 



in Eq. (2.2). For example, for a vanishing magnetic field and for small x (type I region) the transition 
is of first order, while for large x (type II region) the precise characteristics of the continuous transition 
have not yet been completely understood. These fluctuation effects can be discussed perturbatively 
flj] , but must eventually be studied with lattice simulations |T2, E3| . There is no local order parameter 
to distinguish between the phases, but a number of non-local ones, e.g., the vector (photon) mass ||I2] 
and the vortex tension |b| . 



In general, the expansion parameters in the theory of Eq. (2.2) are, dimensionally, of the form (cou- 
pling constant) /(inverse correlation length = "mass"). Choosing the scalar coupling and the scalar 
mass, corresponds to the standard Ginzburg criterion (/symmetric— /broken) x (scalar correlation length)'^ 
> T. In terms of the dimensionless variables y,x, this means 

y < -Wx"^. (4.1) 

Thus, for mean field theory to be valid, one simply has to be sufficiently deep in the broken phase. 
However, this is not the only expansion parameter: choosing instead, e.g., the vector coupling and the 
vector mass, e§/(7r?7iy), one gets a criterion which is more stringent at small x, 

y < -lOx. (4.2) 

In any case, these are only qualitative criteria and the true size of fluctuation effects can only be 
established by numerical means. In many analogous models, the fluctuations change the mean field 
phase structure completely [p4[ . 

For a non- vanishing field B, much new structure appears even on the MF level, as discussed above. 
What happens in the full theory can, in principle, be solved by the numerical methods presented here. 
Many simulations have been carried out with different simplified models (see, e.g., p5| and references 
therein), but it is not clear if the results should agree with the full locally gauge invariant U(I)-fIIiggs 
theory, preserving all relevant degrees of freedom. 



Let us again consider the canonical ensemble (2.10), in which the magnetic field H is fixed and B 



is determined by the system. If the system is in the broken phase a.t H = 0, and H is increased, the 
system eventually changes to the symmetric phase, which is characterized by a zero photon mass. For 
type I, the transition is of first order, and no other phases are believed to exist. 

II 



For type II, the magnetic field penetrates the system at strong enough fields, H > Hd — e^T/2Ti, 
via vortices, and the photon becomes massless. On the mean field level, the result would be a vortex 
lattice. However, it is clear that this picture changes when fluctuations are included. The interaction 
between the vortices becomes weaker when their distance increases and therefore at small enough 
magnetic field, the fluctuations will certainly destroy the lattice ordering. On the other hand, when 
the field is large, the vortices start to overlap and the cannot anymore be treated as interacting 
line-like objects. There is experimental evidence 0] as well as theoretical understanding in other 
models [§, ^, |§1 that the lattice existing at intermediate fields, then melts into a vortex liquid phase, 
which is thought to be smoothly connected to the symmetric phase ||5|]. 

Therefore, it is an important question whether the lattice ordering predicted by the mean field 
theory really exists at all for intermediate magnetic fields in the U(l)-|-Higgs model. If there is no 
lattice phase, then the only true phase transition in the system is at iJci, when the Meissner phase 
changes into a vortex liquid. 

5 Simulations in the type I region, x < 1/2 

The objective of the simulations in the type I region is to verify the qualitative picture discussed in 
Sec. 0, and to determine quantitatively the vortex tension T, the critical field strength He, and the 
surface tension a. The simulations are carried out at one fixed x — 0.0463 and for y = —1,-2,— 4. 
The lattice sizes are A^3 — 16, 28, Ni — N2 — 16, 24, 28, 32. The number of flux units varies from to 
about 33, corresponding to a range of -8/63 from to about 10. The mean field values of the scalar 
and vector correlation lengths are, at (3c — 4, 

:/3Ga = 2a, (y = -2), 



mr V^2^' 



,MF 



— /3Ga«0.6a, (2/ = -2), 

y 



and indicate that the approach to the continuum limit must be checked explicitly (since the vector 
correlation length is only of the order of the lattice spacing). From |2^, we know that the tree- level 
estimates for the masses are reasonably good everywhere in the broken phase. 

The numerical simulations, as presented here, were performed with the use of hybrid Monte Carlo 
update schemes for the gauge and the Higgs fields on parallel computer architectures. These were 
supplemented by simulations on workstation clusters, where the gauge fields were evolved using the 
more conventional 3-hit Metropolis update. A typical simulation for the calculation of the free energy 
change dF{n) consumed (l...several)xlO^ Monte Carlo sweeps of either algorithm. Such values of the 
statistics allow for relative errorbars at a few percent level in the measured quantity. Errorbars have 
been calculated with jack-knife methods and whenever necessary statistical error propagation was 
used. For a more detailed account of the algorithms, error calculation and the calculation of integrals 
we refer to |^ . 

All the measurements are based on measuring the function W{m) in Eq. (^.17). Since all the 



analysis is based on its integrals, i.e. free energies, we show only one example of W{m) itself in the 
type I region in Fig. 0(a); all other cases are very similar. Fig. 0(a) exhibits in a striking way all the 
features associated with increasing the number of flux units as described in Section p. F or smaller to, 
W{m) contains a series of peaks between integer values ra = n. We know from Eq. ( |2.18| ) that at each 
integer W{rn) = 2m/NiN2. This value is reached within errorbars (not shown for clarity). 

The peak in between the integer values represents the addition of one more vortex line, and the 
integral over the peak gives the additional free energy. The clustering of flux lines as a cylinder (a 
bubble on the 2-dimensional transverse plane) is seen as a decrease of peaks (Eq. ( [3.12| )), and the slab 
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n^l W(m) 



(a) 



(b) 



(c) 




nrr. W(m) 




npf W(m) 




Figure 2: The function 'K(3yW{m) (see Eq. ( p.!?!) ) measured for (a) x = 0.0463 (type I), y = -2, 



/3g = 4, on a 18^ x 16 lattice; (b) a; = 2 (type II), y = -4, /3g = 4, on a 32^ x 16 lattice; (c) x = 2 



(type II), y = -1, /3g = 4, on a 16^^ lattice. The straight Hne is W{m) = 2m/{NiN2) (cf. Eq. ( |2.18| )) 
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Figure 3: T and Ta as a function of iVi = N^ dX (3g = ^, x = 0.0463, y = -I. 



as a relatively constant series of peaks (Eq. ( 3.13| )). At some large n the peaks suddenly disappear 
and the system goes over to the symmetric phase. At this point it may be illuminating to look also at 
the corresponding figures for the type II region in Figs. Q(b,c), which do not seem to show any abrupt 
transition. 



5.1 Vortex tension and interaction energy from n = 1, 2 

We begin the study of the integrals of W{m) by considering them up to n = 1, 2. The vortex tension 
T was already studied in detail in H , and a comparison of T and T2 carried out here is a measurement 
of vortex- vortex binding energy, which allows to differentiate between type I and type II regions. We 
choose y = —1,-2 and study finite size and finite lattice spacing effects in T and T2- The values are 
obtained from Eqs. (3^), (3.3). 

In Fig. g we display T and T2 as a function of iVi = N2 a.t j3g — ■^ and y — —\. The data shows no 
finite size effects. Fitting to a constant value one obtains 



T/el 
T2/el 



73.0(13) 
64.2(6) 



(y = -i,/3G = 4), 

(y = -l,/3G = 4). 



(5.1) 



The fact that T > T2 indicates an attractive force between the vortices. 

To study the approach to continuum, we display in Fig. T and T2 as a function of Pq — ae\ on 
a 28'^ lattice at y = — 1. We can use a single lattice size since finite volume effects were observed to 
be small. A dependence on a is clearly seen and linear continuum extrapolations (the straight lines 
in the figure) result in 

T/el = 81.9(37) (y = -1,/?^ = 00) 
T2/el = 68.6(30) (2/ = -1,/3g 



CX) 



(5.2) 



According to Eq. (3.2) the mean field value at this x is T^^ = -42.75yei so that T/T^^ = 1.91(8). 



The magnitude of the scaling corrections is at the level of 10% for T, somewhat smaller for T2^ if the 
continuum results are compared with those at (3g = 4. 

To study the y dependence, we have measured T and T2 also at y = —2. The data again shows no 
apparent finite size effects. The result of a fit to a constant value is 



T/el 
T2/el 



118.9(21), 
106.2(8), 



(y = -2,/3G = 4) 
(y = -2,/3G = 4) 



(5.3) 
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Figure 4: T and T2 as a function of /3q^ on a 28'^ lattice and for x — 0.0463, j/ = — 1. 

corresponding to T/T^^ = 1.39(2). Again, finite Pa effects are expected to increase this number 
slightly. 

Summarising this subsection, one observes clearly a finite tension difference T — T2 caused by the 
attractive force between the vortices. There are also clear deviations from mean field which are largest 
close to the transition point, and actually quite sizable at y = —1. 



5.2 Results for large n 

Data for H{B) (see Eq. ( [2.19| )) is presented in Figs. 0(a-c) for three sets of parameter values, y = — 1 
on a 16 X 32^ lattice, y — —2 on a 18'^ lattice, and y — —4 on a 16'^ lattice. We remind that this data 
contains integration over the successive peaks of W{m) such as those in Fig. ||(a). 

In all the figures, in agreement with Section]^ , one can see a region with a rapid variation (cylinder 
region), in which the data can be fitted to Eq. ( |3.19| ), 



H{B) ^Hc + a- 



L1L2BHC 



(5.4) 



and a constant region (slab region), fitted to Eq. ( 3.20 ). The constant region gives the value of He, 
the rapidly varying region t he va lue of a. After these values are known, one can compute the two 
transition values: Bi in Eq. ( 3.15 ) for the cylinder — > slab transition and B2 in Eq. ( 3.17 ) for the slab 
— > symmetric transition: 

i?i = — , S2 = i/c-2e3\/^. (5.5) 



A^i 



These are plotted in Figs, ^(a-c) as vertical lines. One observes that the flat section is approximately 
located wit hin th e B-interval Bi < B < B2- For still larger lattices there would be a second cylinder 
stage (Eq. ( [3.14 )), but now the condition ( 3.1S| ) is satisfied and the transition takes place directly 
from a slab. 

The horizontal lines in all three figures correspond to fits to the flat sect ion. They determine the 
critical field He, which in the mean field approximation is given by Eq. ( |3.8| ), H^^ = — 3.286ye|. We 
obtain 

He = 5.52(3)ei (= 1.680(10)iJ^^P) {y = -1), 

= 9.87(8)e3 (= 1.502(12)i?MF) (2/ - -2), (5.6) 

= 17.17(9)e3 (= 1.306(07)i/,MF) (y - -4). 
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Figure 5: The magnetic field H{B) (Eq. (|2.19D ) for x = 0.0463 as a function of B (Eq. (^): (a) 
y = —1 on a 16 X 32^ lattice; (b) y = —2 on a 18'^ lattice; (c) y = — 4 on a 16'^ lattice. The flux 
number n can be simply counted by starting from the left. For vertical lines, see Eq. (5.5). 
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The curves in the figures are fits to H{B) for values B < Bi with the form ( 3.19 ) and they 
determine a. At small n one sees deviations from a spherical bubble shape and correspondingly the 
fits use y = —1 data with n > 1, y = —2 data wi th n > 4 and y = —4 data with n > 6. Compared 
with the value of a^^ = 4.61(-y)3/2e| from Eq. (|U), the fit results in 

a = 12.7(3)et (= 2.76(7)a^F) (y - -1), 

= 20.2(12)ei (= 1.55(9)aMF) (y = ~2), (5.7) 

= 51.3(30)e4 (= 1.39(8)ctMF) (y = -4). 



Finally, let us consider the thermodynamical discontinuities discussed in Eqs. ([3.22). We note that 



at finite volumes, AS is directly given by the value (w B2) where the system goes into the symmetric 



phase. However, there are large finite volume effects in this quantity (cf. Eq. (3.17)). In the infinite 
volume limit, we expect rather that A_B w H^^ with the values of H^ as given in Eq. (^.6|). Through 
Eqs. (3.23|), (3.24) and the curve Hc{y) from Eq. (5.6), this determines the latent heat. 



Concluding this section, we observe that we have been able to determine T — T2, verifying that we 
are indeed in the type I region, as well as all the main thermodynamical properties of the first order 
transition from the Meissner phase to the Coulomb phase. Qualitatively, the mean field picture is 
valid, but on a quantitative level we observe a transition which is significantly stronger than at the 
mean field level, particularly for small \y\. 

6 Simulations in the type II region, x > 1/2 

The simulations in the type II region have been carried out at x — 2, y = —0.4, —1, —4, —8. The data 
for y = —0.4, —1 is mainly for the n = 1, 2 sectors of the theory to study vortex- vortex interactions 
and their finite size effects on periodic boxes. The data for y — —4, —8 is for large winding numbers 
and serves the purpose of determining the properties of the theory at a finite field H. 

As in the case of type I, all runs first compute W{m) but the analysis is based on its integrals. We 
thus present only two examples of W{m) in Figs. 0(b,c); the rest are similar and differ characteristically 
from those for type I in Fig. 0(a). 

The analysis then proceeds as follows: 

• To have a rough picture of the spatial structure of a single vortex, we measure it in two ways: 
by explicitly measuring the profile of the vortex and by measuring the vector correlation length, 
which is expected to set the scale of a type II vortex at large distances. 



To have an understanding of vortex- vortex interactions we use finite size scaling of T (Eq. (3.29|)), 



the tension of a single vortex, and of T2 (Eq. (3.30)), the free energy per vortex of a two-vortex 



• 



system. We check that in the infinite volume limit T = T2. 

To confirm the understanding so obtained, we check whether the measured free energies of a 
many-vortex system, n = 1,2, ...,~ 50, can be reproduced. We also attempt to determine the 
phase of the system in the many-vortex region. 

6.1 Vector field mass 

To have information on the physical distance scales of the system, we shall first for n = and in 
the broken phase determine the vector field Ai mass my (the "photon mass" ) . We remind that this 
quantity is a (non-local) order parameter for the transition: it vanishes in the symmetric phase even 
after non-perturbative effects are taken into account |l3]. Since we are in the type II region, the 
scalar mass ms is always heavier than my and is thus subdominant at large distances (for explicit 
measurements of ms , see [O] ) . 
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In r(d) 




Figure 6: Photon correlation functions on a 32^ x 64 lattice at y = — 8 (bottom), —4, —1, —0.4 (top) 
as a function oi d = ax^. The curves correspond to the fit in Eq. (|6.2[) with one massive scale. 



The vector mass is determined as follows. We define the plane averaged correlator 



(6.1) 



choosing the lowest non-zero value fci = 2tt/Ni for fci . We then measure and fit the correlator as 

r(ax3) - (0(0; k,)0*{x3; fci)) = c(e-'^^^ + e-(^3--3)). (5.2) 

The frequency mode lo is related to my via the lattice dispersion relation 



2 2 

TOya 



2 sinh(— oj) 



2sin(-fci) 



(6.3) 



We have measured mv for y = —0.4, —1, —4, —8 and for lattices of geometry A^i = N2 and 
N3 = 2Ni with A^i = 16, 24, 32. In Fig. ||, examples of the correlation functions are displayed. A 
fit in accord with Eq. (6.2) gives a perfect description of the data, and higher states appear to be 
strongly suppressed. 

All values of mv determined here are collected in Table 0. For comparison, the numerical value 
of the mean field mass m^^ — e\^ —yjx is also given in the table. The masses exhibit a mild and 
partly non- monotonous dependence on A^i. This means that we are not yet at large enough volumes, 
particularly for y = —0.4, —1. To have an estimate of the infinite volume values, we use a fit of the 
form 



mv(N-\) — mviNi — 00) + ce 



MF 



Ni 



(6.4) 



assuming finite size effects to be those of a massive theory. Other analogous fit forms give comparable 
results within errorbars. In case of a non-monotonous mv{Ni), only the two largest lattice sizes are 
included in this fit. The boldfaced numbers in Table W correspond to the extrapolations. They will 
be used in the following sections of the paper. 

6.2 Profile of the vortex 

The basic idea of our method of imposing a fiux of B on the system is to put a defect along an arbitrary 
line and let the system dynamically decide how it responds to this, where it places the vortices, etc. In 
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MF 


MF 
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Table 1: Results for the vector mass my. The values in italics are the mean field (MF) results and 
the boldface values are the infinite- volume extrapolations. 

a,Jr)/a,Jr = 0) 




Figure 7: The profile of the flux as a function of the distance r from the center of the vortex. The 
curves {y — —1, top; j/ = — 4 in the middle; y = —8, bottom) correspond to fits according to Eq. (6.7) 
for values of z > 2. The parts of the fit curves corresponding to z > 2 are drawn with solid lines. 



the free energy measurement of a single vortex, its position never enters. How can one then measure 
the profile of a vortex, i.e., the distribution of the flux density ai2 on its planar cross section? 

To do this one has to fix the position of the vortex. This can be done by making the infinitely long 
defect in Eq. ( 2.14 ) of finite length in the X3 direction. This effectively places a monopole and an 
antimonopole at the ends of the defect and localises the response. Concretely, take a 17 x 17 x 32 lattice 
and shift the plaquette values by 012(2:1, 2^2 : 2:3) — > ai2(xi,X2,X3) -I- 27r for fixed values a;i = 3:2 = 8 



but only for values of x^ — 8, 9, ...,23, 24. Then the monopoles are located at 



is 



;-f i,7-|-i) and 



^, 24 + i), which gives a total distance between the monopoles of dmm = 17a = 4.25/e| in 
the a;3-direction. The Dirac string between the monopoles is again compensated for by a vortex, and 
the field profiles are measured at half-way distance between the monopoles, i.e., at X3 = 16. 

Let r be the distance from the center of the defect {8 + ^,8 + ^). Assuming rotational invariance 



19 



we measure the statistical average of the quantity 



ai2ir = 0)' 



(6.5) 



Results are displayed in Fig. M for y := — 8. We have also repeated the procedure for a 33 x 33 x 64 
lattice and the data at y = —4, —1 in Fig. Hstem from simulations with that lattice size. 

In physical units, we expect the vortex size to be related to the vector correlation length 1/my. We 
thus introduce the dimensionless variable 



z — myr, 
and we fit the large distance behavior of the vortex profile for values z > 2 by the form 



ai2{z) 
ai2{z = 0) 



Ci +C2Ko{z). 



(6.6) 



(6.7) 



We find that the asymptotic distance behavior at z > 2 is perfectly described by the fit. Thus the 
asymptotic profile is governed by the vector mass alone. In the cases y = —8, —4 the constant ci can 
be chosen to be exactly zero, while at y = — 1 it has a slightly positive value, which we suspect to be a 
finite size effect. Note that the mean field value for the scalar mass m^^ = ^J—2ye\ dX x — 2 exceeds 
the photon mass by a factor of two, so that its contribution can be neglected. 

Clearly, this analysis is only valid if the distance between the monopoles is larger than the width 
of the vortex, dmm > l/my- The distance should not be too large, either, or the fluctuations of 
the vortex spread the fields and in the limit of infinite separation, the field distribution becomes 
uniform. We estimate that this leads to a Gaussian probability distribution of width ^ y'dn 



the position of the center of the vortex, and therefore, if r > yd^ 
In our case z > 2 guarantees that this is satisfied. 



^/T for 

i/T, this effect can be neglected. 



6.3 Vortex-vortex interactions for n = 1,2 



As discussed in Sec. 3.2, a determination of T, T2 from Eqs. (3.2), (3.3) is subject to much larger 



finite size effects in the type II case than in the type I case. (For an explicit comparison, see Figs. H, 
H) Thus, we have to employ the ansatz discussed in Sec. 3.2 to determine T, T2. We denote the finite 
volume values by T{L),T2{L). 

All of our data for T{L) and T2{L) at x = 2,y = —l,/3c = 4 measured on cubic N^ lattices 
(L = No) is shown in Fig. @ as a function of the lattice size A''. The pattern is similar for other values 
of y. We recall that the tensions (Eqs. (p.2|),(3.3)) are obtained by measuring W{m) (Eq. (|2.17)) and 



integrating (Eq. (2.19)) over one or the two first peaks of functions such as those in Fig. 

The data in Fig. H is then fitted with two free parameters, namely T and X according to the finite 
size scaling laws in Eqs. (3.29) and (3.3C), using my as determined in Sec. 3.1. The numerical value 
of the distance control variable z is, for A^i = N2 



N. 



z — myL 



0.564:N 


{y = 


-8), 


0A21N 


{y = 


-4), 


0.218A^ 


{y = 


-1), 


0.127A^ 


(2/ = 


-0.4) 



(6.8) 



We include data with z > 1 in the lattice sums in Eqs. ( 3.27 ). Typical finite size scaling studies in 
other lattice models use values of z > 2, but we have checked that our fit results are stable under 
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Figure 8: T'2(L) (triangles) and T(i) (circles) at x 
lattice as a function of N. 



2, Pc ~ 4 and y = —1, measured on an N^ 



further omission of small z data. The fits for T (such as the lower curve in Fig. 0) have Xdot ^ ^.4 
at y = — 1 and x^^j = 0.69 at y = —0.4. The parameters are 



T = 



17.87(35)e§ i 
9.331(91)eii 
2.571(39)e|i 



1.077(21)TMF) 
1.125(11)T^*P) 
1.240(19)T'^P) 
= 1.089(30)e§ (= 1.313(36)TMP) 



X -... (y = -8), 

(y = -4), 

= 1.02(6)[10] (y--l), 
= 0.80(3)[20] (y = -0.4), 



(6.9) 



where the mean field result is T^^ = — 2.0735ye|, and the numbers in the square brackets correspond 
to changing the (infinite volume) values of my within the errorbars given in Table w. For y = —8, —4, 
the data is not good enough for determining X, while for all values of T, the errors from the uncertainty 
in mv are within the statistical errorbars. Both T{L) and T2{L) give comparable values for the 
infinite volume extrapolation (see Fig. 0). It is quite impressive that such large finite size corrections 
are described so well by the fit. This is a clear indication of the fact that the interaction is indeed 
described by Kq- 

Consider then the approach to continuum, a = l/{f3ae'^) -^ 0, oi T and T2. We choose y — —1. 
For a fixed N"^ = 28^^ lattice we determine T and T2 at Pa = 1, 2, 4, 6 and 8. Table |2| contains in its 
second and fourth column data for the measured quantities T{L) and T2{L), L = 2 8a. This data is 
then extrapolated to infinite volume u sing t he pre viou s measurements of X in Eq. (|6.9|) and mv in 
Table 1^, and assuming the validity of ( 3.29 ) and ( 3.3C| ). The extrapolations are given in columns 3 
and 5 of Table g. Note that we have assumed that X , ray do not depend significantly on /3g', and 
we use the same values for each j3g\ this clearly introduces some extra /^c-dependence, but it turns 
out that the outcome of this rough procedure is numerically quite satisfactory. The data somewhat 
scatters; at large values of /3g autocorrelations within the Monte Carlo process are large. 

We display T{L = 00) and T2{L = 00) as a function of ae§ = 1//3g in Fig. g. The straight line fit 
(the curve in the figure) corresponds to the form 



T{(3c 



for T and a similar one for To. The fit results in 



-i 



T = 2.594(83)6^ (= 1.251(40)r"^) {y = -1), 



(6.10) 



(6.11) 
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Figure 9: The continuum extrapolation of T/T^^ and T2/T^^ as a function of 1/Pg at a; = 2 and 
y = —1. Here T^^ — 2.0735e|. Circles and triangles correspond to infinite volume extrapolations of 
rp jrpMF g^^^ T2/T^^ . The solid circle and the fitted line correspond to the continuum extrapolation. 



Pa' 


nL)/el 


r(L = oo)/ei 


T2{L)/el 


UL = ^)/el 


1.0000 


2.488(170) 


2.488(170) 


2.489(125) 


2.489(124) 


.5000 


2.385(179) 


2.385(179) 


2.457(131) 


2.457(131) 


.2500 


2.485(182) 


2.473(181) 


2.838(112) 


2.735(107) 


.1666 


2.428(232) 


2.312(221) 


3.248(198) 


2.656(162) 


.1250 


2.738(276) 


2.326(237) 


4.187(256) 


2.605(159) 



Table 2: The infinite volume (L = cx)) extrapolations of T/e\ and Tije^ as a function of /3g,^ at x = 2 
and y = —1; L = 28a (see text). For (3^^ — 1.0, 0.5, the dependence on L is within the errorbars. 



having Xdoi — 0.83. This value is completely compatible with the one determined in Eq. (6J) at 



y = —1 and thus we conclude that finite-a corrections are small and apparently under control for the 
values of (3g considered. Fig. g can be compared with that for the type I case. Fig. ^. 

To summarize this subsection, we have verified that at a; = 2, T — r2 < at finite volumes but 
—f for L — > cx). Hence we are in the type H region. We then showed that the infinite volume 
extrapolations T{L — oo) can in turn be extrapolated to a finite physical value in the continuum 
limit. 

6.4 The vortex system at large n 



Consider then the case of large n. As in the type I region, we plot H{B) from Eq. (|2.19| ). Our main 
data is shown in Fig. |l^, which contains the integrals over the peaks in W{m), shown in Fig. || 

A first-order phase transition means that B behaves discontinuously as a function of H . In the 
microcanonical picture, it would correspond to a region in which H is constant, or more precisely, 
would exhibit similar finite-size effects as those discussed in Sec. 3.1. The only region where this 
might be the case, is at small B (cf. Fig. |o|(b)). If that turns out to be true, there is a first-order 
phase transition from i? = or a small non-zero value, to some finite larger value of B. In principle 
the transition could be one from a vortex liquid phase at such small B that the distance between the 
vortices is large enough that the interaction cannot preserve the lattice structure, to the lattice phase; 
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or, more likely in a finite volume, from a vortex lattice phase enforced by the periodic boundary 
conditions, to a vortex liquid phase where the vortices are so close to each other that the average 
fluctuations become larger than the average distance. 

In order to check whether a lattice or a liquid phase is more likely, the data for H{B) is fitted to the 
ansatz in Eq. ( 3.31 ) with both a triangular (T) and a square (S) lattice. We remind that now the ansatz 
is for H ~ dF/{VdB) and that the distance control variable z quantitatively is (Li — L2 — L = Na) 
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(6.12) 



^B/e 

where the numbers in the parentheses correspond to the lattice sizes used in Fig. |l^. The fit is done 
for y = -4, -8 with a fixed interval Bq < B < H^^ , where Boiy = -4) = O.Se^, Ba{y = -8) = 1.6ei 
and H^^ — —ye^. We consider T and X as free parameters in the fit. The results for the parameters 
are given in Table ^ and the fits are shown in Figs. [lO|(a,b). 
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Table 3: Fit results for X and T from the fits in Figs. [l^(a,b), with 0.8ei <B < H^^ {y = -4) and 
1.663 < B < H^ {y — —8), assuming the existence of a triangular (T) or a square (S) vortex lattice. 
The mean field value in the dilute limit is X = 0.76. 

To summarize the results for the parameters at different y- values: for the vortex-vortex interaction 
parameter X{y) we have X{—S) — 0.80(5), X{—A) — 0.90(5). There is only mild y-dependence and 
the results are consistent with those at y = — 1,— .4 in Eq. (|6.S|). The mean field value (applicable 
in the dilute limit) is X = 0.76, independent of y. The tension values T{y) of the fit at large B 
are T{—A) = 10.7(3)e§, T(— 8) — 19.9(5)e|. These values differ from the infinite volume tension as 
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Figure 10: The field strength H{B) as a function oi B a.t x = 2.0, /3g = 4 and (a) y = —8 on a 



16 lattice; (b) y = —4 on a 32 x 16 lattice; (c) y = — 1 on a 16 lattice. The vertical lines in 



20" 

(a),(b) indicate the fit intervals as denoted by Bq and H^^ for the triangular lattice fit in Eq. (3.31) 



As discussed in the text, the good agreement of the triangular fit with the data is not enough to 
guarantee that there really is a lattice structure in the system. 
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Figure 11: Finite size scaling analysis of M{B) = B - H{B) at a; = 2,y = -4 and B/e| ^ 6.2 831 
on A''^ boxes (compare with Fig. |l^(b)). The fitted curve corresponds to the ansatz in Eq. (6.14). It 
extrapolates to a non-zero value of M{B) here. 



determined in the dilute limit in Eq. (6.9) by about 10 percent, which may well be a finite volume 
effect. Finally, the S and T lattice fits coincide within errorbars. 

The fact that the S and T lattice fits work equally well and give consistent results, means that we 
are not sensitive to the average spatial distribution of the vortices. Thus, a vortex liquid phase cannot 
be excluded, either. 

To get more information on the phase, let us then inspect the region of large B. At y — —4 and 
y = —8, the peaks in W{m) persist at all values of n studied (see Fig. ||), and the field strength 
H{B) does not reach the mean field Coulomb phase value H(B) = B, corresponding to vanishing 
magnetization M{B) = B — H{B) — 0. On the mean field level the transition to M{B) ~ takes 



place according to Eq. (3.33) at 



B^H. 



MF 



-2/ei, 



(6.13) 

which is shown in Figs. |o|(a-c) as a solid line and is well within the studied range. However, we do 
not observe any structure around these values. 

For the data set sA y — —4, we have performed a finite size scaling study of M{B) at a large value of 
B/el = 6.2831. Considering hypercubic N = Ni= N2=N3 boxes for values iV = 8, 12, 16, 20, 24, we 
display data for M{B) as a function of N in Fig. O. The data differs from the mean field prediction 
M{B) = and suggests an infinite volume extrapolation to a value M{B /e\ = 6.2831) ~ —0.12. The 
fitted curve in the figure corresponds to a finite size scaling ansatz 



M = M{B) 



A 



(6.14) 



with 7=1. Such power law finite size corrections are expected in a massless phase. The exponent 7 
could in principle be non-trivial, but we are currently not in the position to determine its value from 
the fit. In any case, it is clear that the deviation of the measured M{B) from the mean field value 
M = is not a finite size effect. 

When B is further increased, H{B) should finally approach the value B, but achieving high enough 
n is costly in computer time. The required value of n can be reduced by, say, decreasing —y. As an 
example, the y = —1 case is shown in Figs. 0(c) and ^Q(c). Here H{B) reaches B within errorbars 
already at n w 3 but the data is not good for any quantitative analysis. However, it illustrates 
an interesting point: contrary to what one might conclude on the basis of the plots for W{m) in 
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Figure 12: The phase diagrams for x = 0.0463 (left) and x = 2.0 (right). The critical values of He in 
the type I region are from Eq. (p!q ), for Hd in the type II region from Eq. (6.9) (cf. Eq. (3.32)), and 
for y at the point H = from [[12[. 



Figs. g(a,b), on finite lattices there are no cusps at integer values of m; in addition, W{m) can fall 
below the line 2m/{NiN2). In fact, M ^ B — H{B) « is obtained when positive and negative 
contributions cancel out. 

In experiments with high-temperature superconductors 0] a weakly first-order phase transition has 
been found at a large magnetic field H, and it has been interpreted as a melting transition of the 
vortex lattice. Like the transition we could imagine to observe at small B, it would also show up in 
our data as a constant part in Fig. ^ Clearly, there is no sign of such a transition in our data, so 
if one exists in the U(l)-fIIiggs model at the value x = 2 we have studied, the transition must be 
extremely weak. A much more natural conclusion seems to be that at least for the parameter values 
studied, the 3d U(l)+IIiggs model does not have a vortex lattice phase as a thermodynamical state, 
but that the structure obtained at the mean field level is removed by thermal fluctuations. Thus we 
would be in the liquid phase all the time, and it just smoothly changes into the symmetric Coulomb 
phase. This is in agreement with the behaviour of my as a non-local gauge invariant order parameter: 
there is only one "Meissner" phase at H < Hd = T/(2tt) with my ^ 0, while for H > Hd, my 
vanishes. 



7 Conclusions and outlook 

In this paper, we have shown how the interactions of vortices can be studied non-perturbatively with 
lattice simulations, starting directly from field theory. We have shown that the interaction energy 
between two vortices is negative in the type I region and positive in the type II region, as expected 
from mean field (MF) theory. Moreover, increasing the magnetic field further on, we have shown that 
vortex interactions can lead to phases with different types of macroscopic structures. 

The present simulations were carried out in a region i n wh ich one approaches the limits of validity 
of MF-theory. In superconductivity —y ^ 1 (see Eq. (2^)) and MF-theory works very well. Our 
simulations span the range y — —8, ..., —0.4 and we, indeed, find increasing deviations from MF. 

The regions of the phase diagram displaying Meissner, vortex liquid, vortex lattice, and normal 
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Coulomb behaviour, have to date not been accurately mapped out in the full Ginzburg-Landau theory. 
We have demonstrated that the methods introduced here can answer these questions. Our present 
results and conclusions for the phase diagram in the (?/, iJ)-plane at x = 0.0463, x = 2 are shown in 
Fig. |l^. At a; = 0.0463, we find a first order phase transition in qualitative accordance with the MF 
estimate. However, the transition is much stronger than at the MF level, particularly for small He 
(in fact, at the MF level the transition is of the 2nd order for He -^ 0, while the actual transition 
continues to be of the first order there). The vector field mass is non-zero only in the Meissner phase. 

In the type II case x — 2 and at y = —4, we observe only a single transition from the Meissner phase 
to a vortex liquid phase, which then smoothly turns into the Coulomb phase. This is in strong contrast 
to the MF prediction: we do not find any indications for a phase transition associated with the melting 
of a vortex lattice phase at large B « H^^ , so it seems natural to conclude that we are in the liquid 
phase all the time, and the fluctuations are strong enough to remove the lattice structure predicted at 
the MF level. On the other hand, wc found very preliminary indications (Fig. |l^(b)) that the transition 
to the vortex liquid phase at small B might be of the first order. A possible interpretation for this 
behaviour is a finite volume effect associated with the 'melting' of an unphysical lattice structure 
induced by the periodic boundary conditions, for a small number of flux quanta (in other words, a 
transition from a boundary-effect-dominated region to a bulk-dominated region). 

Whether or not the theory in the thermodynamical limit allows for a physical triangular vortex 
lattice state in a strict sense at large values of — y, currently is an open question. It remains an 
algorithmic and computational challenge to extend the present simulations to the case of a truly 
macroscopic number of vortices. 

One could imagine that techniques analogous to those introduced here can be developed for defects 
other than vortices, as well. For instance, in a theory with monopoles, one could flx the total magnetic 
flux emerging from the whole lattice volume. This would allow for a gauge-invariant non-perturbative 
study of the free energy of one monopole or of an ensemble of several monopoles. 

Finally, let us note that the present techniques can be easily extended to the 3d SU(2)xU(l)-|-Higgs 
theory, relevant for the cosmological electroweak phase transition in the Standard Model and many 
extensions thereof. In that case, an external magnetic field can represent physical (cosmological) 
reality, and can in principle lead to the emergence of phases analogous to the vortex lattice of type II 
superconductors |2^ . Such phases might be relevant for baryogenesis and cosmology, and would thus 
be of considerable interest. We have observed that even in the U(l)+IIiggs system where the lattice 
prediction is much more robust to begin with, fluctuations tend to remove the structure and result 
in a vortex liquid phase. This makes it quite understandable that in the SU(2)xU(l) case no lattice 
structure (nor, in fact, a clear liquid-like behaviour) could be observed in the region of the parameter 
space studied so far B, but one rather finds a symmetric phase. 
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